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1 Introduction 

The virtual testing of delamination is a goal shared by many practitioners, especially in the aeronautical 
field. In order to reach such an objective, two research topics which have undergone drastic changes 
over the last twenty years must be linked: the relevant modeling of composites and the efficient analysis 
of structures. 

Indeed, there have been many advances toward a better understanding of the mechanics of lami- 
nated composites and the mechanisms of damage. The validity of two types of models, microscale mod- 
els and mesoscale models, has been proven. Microscale models are closely connected to the physics of 
the material and, thus, provide a reliable framework for simulation. They take into account many dam- 
age processes, such as diffuse intralaminar degradations percolating into transverse cracking, diffuse 
interface degradations leading to distributed delamination, chemically- or thermally-induced degrada- 
tions, or fiber breakage. |12l IB"]. On the microscale, simulations can combine continuous (damage) 
and discrete (fracture) degradation models [18 . Unfortunately, the analysis of models defined on 
the microscale requires such a refined discretization that only small test specimens can be simulated. 
Industrial-size structural calculations are beyond the reach of even recent computers. Mesomodels 
[TJ H31 E] are defined on a scale which makes both the introduction of physics-based components and 
the simulation of small industrial structures possible. Very often these models rely on the definition 
of two mesoconstituents, the ply (a three-dimensional entity) and the interface (a two-dimensional 
entity), which are modeled using continuum (damage) mechanics and behavior derived from the ho- 
mogenization of micromodels |171 116) . Nevertheless, in order to achieve reliable simulations, refined 
discretizations are still required for the correct representation of the stress gradients induced by edge 
effects, which are responsible for the initiation of many degradations. Therefore, the resulting problems 
remain very large (in terms of the number of degrees of freedom) and highly nonlinear, which creates 
potential instabilities. 

In a first approach to the reliable simulation of delamination in composite structures, we chose 
to neglect the effect of damage within plies and concentrate on the degradations at the interfaces. 
Thus, we adopted the mesomodel presented in [T], in which the debonding phenomenon is localized at 
the interfaces and handled through cohesive behavior. A similar approach with a different interface 
behavior (degradations based on plasticity) was applied in |22) . 

In order to handle the large nonlinear systems associated with this modeling approach, one can 
consider using one of the several multiscale [HJ |5J 0J [T3] and enrichment [2"Tl ITUI |H1 HO] techniques 
developed recently. We based our strategy on the mixed domain decomposition method described 
in |14j . which places special emphasis on the interfaces between substructures. Consequently, the 
reference problem resulting from the mesomodel chosen is substructured by nature, and the cohesive 
interfaces of the model are handled within the interfaces of the domain decomposition method. This 
idea is developed in Section [2] Furthermore, the resolution of the substructured problem by a LATIN 



iterative solver has very advantageous numerical properties: the nonlinearities are dealt with through 
local problems, and very little matrix reassembling is required. The incremental micro-macro LATIN 



algorithm as a resolution strategy for delamination problems is presented in Section 2.2 As shown in 
Section 2.3 the direct application of this method leads to a number of numerical difficulties. A first 
issue occurs when setting the parameters of the method: in Section [3j we present the indispensable 
tuning of the search directions according to the interface's status. Subsequently, the main remaining 
difficulties concern the treatment of the macroscale of the problem. In this paper, the emphasis is 
on the adaptation of our strategy in order to deal with large macroproblems. (Important remarks 
on how to make the macroproblem more relevant can be found in In Section |4j we present the 

parallelization of the resolution, which was inspired by [19j . In order to do that, we introduce a third 
level of discretization after the first and most refined level (the finite element) and the second level (the 
substructure) : interconnected substructures are combined into "super-substructures" (which fill up the 
memory capacity of processors) connected to one another through "super-interfaces" using Message 
Passing Interface (MPI). The method is validated in Section [5] using a complex test case. The handling 
of the instabilities is not described in this paper. The interested reader may refer to [11] where the 
adaptation of an arc-length algorithm with local control is presented. 



2 Application of the two-scale domain decomposition strategy 
to delamination analysis 

2.1 The substructured delamination problem 

Let us consider a laminated structure E defined in a domain fl bounded by d£l and consisting of 
Np adjacent plies P, each defined in a domain fip, with ft = [J p ilp. Adjacent plies P and P' are 
joined by cohesive interfaces Ipp>. An external traction field F_ d is prescribed over a part <9f2/ of f2, 
and a displacement field U_ d is prescribed over the complementary part dQ u of f2. Let n P denote the 
outer normal to the boundary d£lp of Ply P, / the volume force, a the Cauchy stress tensor and e the 
symmetric part of the displacement gradient. The simulation is performed using a classical incremental 
scheme, assuming small perturbations and quasi-static isothermal evolution over time. 
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Figure 1: Decomposition of the laminated composite structure into substructures 
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The laminated structure E is decomposed into substructures and interfaces as shown in Fig. [T] 
Each of these mechanical entities has its own kinematic and static unknown fields as well as its own 
constitutive law. The substructuring pattern is defined in such a way that the domain decomposition 
interfaces coincide with the material's cohesive interfaces, so that each substructure belongs to a unique 
ply P and has a constant linear constitutive law. A substructure E defined in Domain VLe is connected 
to an adjacent substructure E' through an interface Tee' = dflE H 80,e' (Fig. [2]) . The surface entity 
Y ee' applies force distributions F E , Fe 1 an d displacement distributions W E , W E > to E and E' . Let 
Te = Ub'geTbb'- Over a substructure E such that Te D <9f2 ^ 0, the boundary condition (£7 d ,F d ) 
is applied through a boundary interface T E d ■ 



Ue,cf e , 



(F e , W e ) 
(F E >,W E > 




Figure 2: Decomposition of the laminated composite structure into substructures: mixed unknown 
fields 

At each step of the incremental time resolution algorithm, the substructured quasi-static problem 
consists in finding s = (se)eee (where se — (Mb; W e , & e < Fe)) which is a solution of the following 
equations: 

• kinematic admissibility of Substructure E: 

at each point of Te, Ue = W. e (1) 

• static admissibility of Substructure E: 

V(u E *,W E *) eUexWe I u E * ]aQB = W E \ 
J n Tr {g E ^E*i) d n = jf l d .u E *cM (2) 



F E .W E * dT 



linear orthotropic constitutive law of Substructure E: 



at each point of £Ie, Q_ e = Ke(« £ ) 



• behavior of the interfaces Tee 1 € Te- 



at each point of Te e' € Te, 
Ree>(W e ,W e „F e ,F e ,)=0 



(3) 



(4) 
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behavior of the interfaces at the boundary dfl nTg: 



at each point of T Ed , Rf.AW E , F E ) = 
(W E = u d on d£l u and F E = F d on <9f2j ) 

We make the formal relation R EE ' =0 explicit in the two cases we will be considering: 
• perfect interface: 

' F R + F R , =0 



(5) 



cohesive interface: 



F E =K p {[W] EEI {r<t))\W] EEI 
F E + F E ,=0 



(6) 



(7) 



where [W] = W E , — W E . The stiffness operator K_ p can be expressed in the (N_i, N_ 2 ,n P ) 



basis as (Fig. 



/ (l-di)fcl 









(1 - d 2 )fc 2 ° 

(l - h+[W] EE , ds) fc 3 ° y 

where /i+ is the positive indicator function. 



V 



(8) 




Figure 3: The components of the mesomodel 



The cohesive constitutive law of an interface joining two adjacent plies is described classically 
through continuum damage mechanics. Local damage variables (<ii)igri 3], with values ranging from 
(healthy interface) to 1 (completely damaged interface), are introduced into the interface model in 
order to simulate its progressive softening. The parameters (<i;)je[i 3] are related to the local energy 
release rates (5i)jgfi 3] of the interface's degradation modes (traction along Direction N_ 3 and shear 
along Directions jV^ and N_ 2 )- 



Yi 



de d 
ddi 



where 



Y 2 
Y 3 



2*2 (M 



2 6 



EE' 
EE' 



.N 



(9) 



The damage variables are assumed to be functions of a single quantity: the maximum over time Yu of 
a combination of the energy release rates (^It),^ 3 i T<t : 



Y\ t = sup [T < t) (y 3 ? t + 7l Y 1 f T + l2 Y 2 f T 



(10) 
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The evolution laws express that: 



n [Y\ n 

d\ = d 2 = e? 3 = w(Y) where, in general, w(Y) = I — ) (11) 

n and a being scalar parameters of the model. When Parameters 71 and 72 are set to identified 
physical values such that 71 ^ 72 7^ 1, the energies dissipated during the propagation of the crack are 
different for the three modes. Details on the identification procedure for such a model can be found 
in0. 

After a cohesive interface has become fully damaged, it is converted into a (frictionless) contact 
interface. 



2.2 Two-scale iterative resolution of the substructured problem 
2.2.1 Introduction of the macroscopic scale 

In the end, the substructured problem defined in the previous section will be solved using an iterative 
LATIN algorithm, which will be described in the next section. In order to ensure the scalability of the 
strategy, a coarse global problem, associated with the equilibrium and continuity of what one calls the 
"macro" force and displacement fields at the interfaces, must be solved at each iteration. 

Over each interface T EE i such that (E,E') s E 2 , the interface fields are divided into a macro part 
(superscript M) and a micro part (superscript m). The macro part belongs to a small subspace (9 
macro degrees of freedom per plane interface for a 3D problem) . 



F E = F E ' + F% 

w E = w¥ + w r i 



(12) 



The macro and micro data are uncoupled with respect to the interface's virtual work: 

V(F E ,W E ) eT E xW E , 

F E .W E dT= [ F^.We 1 dT+ f F^.W^dT ^ 

Each macrospace is defined by one's choice of its basis. Numerical tests have shown that the use of a 
linear macro basis gives the method good scalability properties. Indeed, the corresponding macrospace 
includes the part of the interface fields with the largest wavelength. Consequently, according to Saint- 
Venant's principle, the micro complement resulting from the iterative resolution of the local problems 
has only a local influence. 

2.2.2 The iterative algorithm 

Here, the iterative LATIN algorithm for the resolution of nonlinear problems is applied to the resolution 
of the substructured reference problem, the nonlincaritics being localized in the (cohesive) interfaces. 

The equations of the problem can be divided into two groups: 

• linear equations in substructure variables and interface macroscopic variables: 

— static admissibility of the substructures 

— kinematic admissibility of the substructures 

— linear constitutive law of the substructures 

— linear equilibrium of the macro interface forces 

• local equations in interface variables: 
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— behavior of the interfaces 

The solutions s = (se)e£E = (u E , W E , a , F E ) r p f, of the first set of equations belong to Space 

Ad and the solutions 's = (se)sge = W E , o , F E ) bpf, of the second set of equations belong to 
r. The converged solution s re f is such that: 

S re /GA d p|r (14) 




Figure 4: Illustration of the LATIN iterative algorithm 

The resolution process consists in seeking the solution s re f alternatively in these two spaces: first, 
a solution s n is found in Ad, then a solution s n+ i is found in T. In order for the two problems to be 
well-posed, one introduces two search directions, E + and E - , linking the solutions s and J through 
the iterative process (see Fig. [4]). Hence, an iteration of the LATIN algorithm consists of two stages: 

• a local stage: 

Find's n+ i € r such that (^s n +± ~ € E + (15) 

• and a linear stage: 

Find s n+ i £ Ad such that ^s„ +1 — ^n+|) "= 0-®) 
In the following sections, the subscript n will be omitted. 



The local stage During the local stage, uncoupled problems are solved at each point of the interfaces 
(Tee')\(e,e')ge 2 ( as wen as (r_E d ),EeE for the interfaces which belong to the boundary cMl): 



Find (F E , W_ E , F_ E , , W E , ) such that: 



■R EE ,(F E ,W E ,F E ,,W E ,)=0 



F, 



(F E ,-F E ,) 



k (Wjs - 
~k+(W E , 



We) 
- W T 



(17) 







The last two equations of this system define the search direction E + (k + and k are scalar search 
direction which, physically, are analogous to "stiffnesses"). In the case of a cohesive interface, Problem 



(17) is nonlinear and its solution is obtained through a Newton- Raphson scheme. 



G 



The linear stage The linear stage consists in the resolution of a series of linear systems within the 
substructures under the constraint of macroscopic equilibrium of the interface forces. 



at Interface T EE '\{e,e')^^ F A J + F% = (18) 



In order to verify the macroscopic condition exactly and the search direction E defined in ( 16 ) as 
well as possible, we use a Lagrangian formulation, which leads to: 



vr £ ew & J (f e -f e ).w e * dv 



J (w E - W.e) - k~w M ^j .w E * dr = o 



(19) 



~M 

which can be viewed as a modified search direction, the Lagrange multiplier W_ becoming an addi- 
tional unknown of Interface T EE > . 



The problem which needs to be solved for each substructure E is obtained by substituting (191 
into (2): 

y{u* E ,W* E ) GU E x We 



j Tr(g(u E ) Kg(u E *)) dil + I k-W E .W E *dT 

f..u E * dn + {F E + k-W E + k~W ).W E dT 



(20) 



n E - d Jr 



The condensation of this equation onto the macro degrees of freedom leads to a relation between F B 

-Af 



and W_ E which can be introduced into the macro equilibrium equation (18 1. Finally, one gets a small 
linear system defined in the macro degrees of freedom. All the subdomains contribute to that "global" 
system through homogenized (condensed) flexibilities h E calculated explicitly. 



~M* x ^ f — M ~M* 

V W £ W M , ^2 h E W .W dT 

E 

= / Fd-W dT-J2 Fe-W dT 



(21) 



The macroscopic problem is discrete by nature and is expressed in matrix form as L M W M — F AI , 
where W M is the vector of the components of the Lagrange multiplier in the macro basis. 



The right-hand side of Equation (21 1 can be viewed as a macroscopic static residual from the 



calculation of a single-scale linear stage. In order to derive this term, Problem (20) must be solved 

independently within each substructure. The resolution of the macroscopic problem (21 1 leads globally 

~M 

to the Lagrange multiplier W_ , which is finally used as a prescribed displacement for the resolution 



of the substructure- independent problems ( 20 1 



In order to carry out the resolutions of ( 20 ) in substructure variables, one uses the finite element 
method. Since the constitutive law of the substructures is linear, the stiffness operator of each sub- 
structure can be factorized once at the beginning of the calculation and reused without modifications 
throughout the analysis, which makes the method numerically advantageous. 

Algorithm [T] summarizes the iterative procedure described in this section. 
2.3 First example of a delamination analysis 

A first example of quasi-static delamination analysis is shown in Fig. [5j The structure is a [0 90] s 
double cantilever beam (DCB). The loading leading to Mode-I quasi-static propagation of the crack is 
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Algorithm 1: The two-scale domain decomposition solver 



Construction of each substructure's operators ; 

Calculation of each substructure's macro homogenized behavior Ljjj ; 
Global assembly of the macroscopic operator; 
Initialization Sq £ T; 
for n = 0, . . . , N do 

Linear stage: calculation of s n G Ad ; 

□ Calculation of the macroscopic right-hand term for each substructure ; 

□ Global assembly of the macroscopic right-hand term ; 

□ Resolution of the macro problem ; 

□ Resolution of the micro problems ; 
Local stage: calculation of s n+ i e T ; 

□ Resolution of the local problems at the boundary interfaces (T c i)e£E 

□ Resolution of the local problems at the interfaces (Tee')(e,e')ge 2 ', 
Calculation of a global error indicator 

end 



increased linearly through ten time steps, the first two corresponding to the initiation of the delami- 
nation and the remainder to the crack's propagation. 

The calculations were performed using a CH — h implementation of the mixed domain decomposi- 
tion method capable of handling the quasi-static analysis of 3D nonlinear problems. In this code, the 
parallel computations use the MPI library to exchange data among several processors. Each proces- 
sor is assigned a set of connected substructures and their interfaces; then it calculates the associated 
operators and solves the local problems (Fig. [8]). The allocation of the substructures among the CPUs 
is handled by a METIS routine. The resolution of the macroproblem does not take full advantage of 
the parallelism because the substructures send their contributions to the macro problem to a separate 
processor in which the matrix is assembled and factorized and the substitutions are performed. 

The direct use of the multiscale domain decomposition strategy to simulate the DCB case led to a 
number of numerical difficulties: 

• The convergence rate of the LATIN-based strategy is highly dependent on the residual stiffnesses 
of the cohesive interfaces as well as the values of the search direction parameters. The iterative 
solver is even likely to stall when using low values of the search direction parameters. In the next 
section, we will briefly describe a practical tuning algorithm for the strategy which guarantees 
convergence. 

• The method looses its numerical scalability when the crack's tip propagates. This phenomenon 
appears clearly in Fig. [6] under the label "No subresolution" . When the delamination process 
propagates (time steps 3 to 10), the number of LATIN iterations required for convergence becomes 
very large. A solution to this problem, described in [IT], enabled us to recover the scalability 
of the method for our test case (under the label "Subresolutions" ) . This approach, which still 
needs to be generalized, is based on a filtering of the long-range effects of the crack's tip achieved 
through the resolution, at each global LATIN iteration, of local nonlinear problems in a box 
surrounding the front (where the main sources of nonlinearities are located) . 

• In this case, the ratio of the number of microscopic DOFs to the number of macroscopic DOFs is 
relatively small (40) . The direct resolution of the macroscopic problem would become an issue if 
one were addressing the simulation of a realistic composite structure. A solution to this problem 
is discussed in Section |4j 
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Figure 5: The four-ply DCB test example 



3 Analysis of the parameters of the iterative algorithm 

A necessary condition for the algorithm to converge is for the search direction parameters k + and k~ to 
be positive definite, symmetrical operators. Previous studies have shown that there exists an optimum 
set of these operators. However, the optimum values are known to be difficult to interpret when the 
interface constitutive laws are complex, and even in simplified cases (perfect interfaces) are expensive 
to calculate. Therefore, our objective was to derive an efficient scalar approximation of these search 
direction operators for debonding analysis purposes. 

As explained in [11 J , the non-monotonic relation between the interface stresses and the displacement 
gap due to damage imposes restrictions on the choice of Parameter k + . Concerning k~ , optimum values 
for a given damage map can be found after a micro/macro decomposition of this operator. Since the 
status of the interfaces changes with the evolution of the delamination (from elastic to damaged, then 
from damaged to ruined) , the search direction parameters need to be updated often in order to remain 
optimum. Therefore, parameters whose efficiency range when they are not optimum is broad enough 
to require less frequent updating are preferred. An effective practical choice is the following: 

• Parameter k + : in order to avoid stalling or divergence of the algorithm, this parameter is set to 
a very high value (i.e. the search direction E + is quasi-infinitely stiff). 

• Parameter k~: 

— perfect interfaces: k~ is set to the classically recommended value E/L [14] . where E is the 
Young's modulus of the adjacent substructure and L a characteristic length of the interface. 

— interfaces with prescribed forces (respectively displacements): k~ is set to a very small 
(respectively large) value in order to enforce the boundary condition through penalization 
in the adjacent substructure. 

— cohesive interfaces: k~ is set to the stiffness of the undamaged interface. 

— delaminated interfaces: k~ is set to zero in the shear direction. In the normal direction, fc~ 
is set to zero in traction and to the initial stiffness in compression. Therefore, the status of 
the interface must be checked regularly (e.g. every ten iterations). 

The use of an infinitely stiff search direction E + and of the initial cohesive interface stiffness as the 
search direction parameter E~ brings us back to a well-known situation. The algorithm can be viewed 
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Figure 6: Subiterations near the crack's tip 



as a secant Newton algorithm in which the solutions of the prediction steps are in equilibrium only in 
the macroscopic space, the equilibrium of the microscopic quantities being achieved at convergence. 



4 The three-scale domain decomposition strategy 

The decomposition into substructures described in Section[2]leads to a very large macro problem and an 
unnecessarily refined macroscopic solution. In order to solve large problems such as that represented 
in Fig. [7J one must place the emphasis on the parallel resolution of the macroproblem and on the 
selection and transmission of the large-wavelength part of the macroscopic solution. 

These two features can be introduced into the method by using any Schur-complement-based do- 
main decomposition technique [TJ. We chose to solve the macroproblem using the BDD method P~^l ITS] . 



4.1 Resolution of the macroproblem through the balancing domain decom- 
position method 

4.1.1 Partitioning of the macroproblem 

The substructures of the initial partitioned problem are grouped into super-substructures {E) separated 
by super-interfaces r j^e' (Fig- [§])• The algebraic problem to be solved within each of these super- 
substructures (dropping Superscript M) is: 



(E) t (E) 
(E) T (E) 



J 66 



w, 



(E) 
(E) 



w h 



(22) 



E 



A (E)~ X (E) 
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(a) Stresses and displacements in the substructures 




where Subscripts b and i refer respectively to the super-interface quantities and to the internal quanti- 
ties of the super-substructures. is a Boolean operator which localizes data in such a way that the 



second equation of System ( 22 1 expresses the continuity of the kinematic unknowns (deduced from a 
single unknown W_ b ), while the third equation expresses the equilibrium of the nodal reactions at the 
super- interfaces. 

First, the local equilibrium is condensed onto the super-interfaces by introducing the Schur com- 
plement and the condensed force F^F ' . The assembled condensed problem becomes: 



where 



SW b = F c 



F c Jj2 A{m £- { c E) 



S^>=lif>-L[f> Lif>" 



F^ = F, 



hi, 

(E) 



■ (B) T (Ey 
J bi ^ii 



■ (E) 



Am 



(23) 



This substructuring technique can be used exactly as in [9] in order to bind a domain which is 
prone to localization and damage to an undamaged region. Nevertheless, for large interface problems 
such as those encountered in our case, the condensed problem is much too large to be solved directly, 
and iterative solvers must be used. 



4.1.2 Resolution of the super-interface problem 

The condensed macroproblem is solved iteratively using a conjugate gradient algorithm. Classically, 
this resolution involves only matrix-vector products and dot products, which are compatible with 
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Figure 8: Three-level substructuring: assignment of substructures to processors 



parallel computation. The recommended Neumann-Neumann preconditioncr S 1 involves the use of 
the pseudo-inverses of the Schur complements of the super-substructures: 

S- 1 =£A^S<^ + A^ T (24) 

E 

The use of this preconditioner means that the inverse of the global super-macro operator is approxi- 
mated by the assembly of the inverses of the local Schur complements. Let us note that the description 
chosen for the interface macrofields precludes the existence of degrees of freedom belonging to more 
than two substructures; consequently, no scaling is required in the preconditioner (at least as long as 
the interfaces are not excessively heterogeneous). The use of pseudo-inverses is associated with an 
optimality condition which ensures that rigid body motions are not solicited (self-equilibrium of the 
floating super-substructures). This condition is verified thanks to a projector which makes the residual 
orthogonal to the kernels of the super-substructures (and possibly to other given subspaces) at each 
iteration of the conjugate gradient. 



4.2 Results 

Fig. [9] shows the convergence rate of the LATIN algorithm when the conjugate gradient algorithm 
for the condensed macroproblem is stopped after a fixed number of iterations. The test case was the 
perforated plate under traction represented in Fig. [7] with the decomposition into super-substructures 
of Fig. [H] It is clear that a rough approximation of the Lagrange multiplier, obtained after very few 
iterations of the conjugate gradient, is sufficient to reach the convergence rate of the multiscale LATIN 
algorithm. Typically, the algorithm is stopped when the residual error (normalized by the initial error) 
falls below 10 _1 . Thus, the third-level enforcement of the admissibility of the macroforces (through the 
projection) appears to be sufficient for the determination of the large-wavelength part of the solution 
to be transmitted through the structure at each iteration of the resolution. 



5 Efficiency of the strategy: study of a complex test case 

In this section, we illustrate the efficiency of the three-scale domain decomposition strategy through the 
simulation of the evolution of debonding in the bolted composite joint shown in Fig. |10| Each composite 
plate interacts with the adjacent plates and with the two steel bolts through contact interfaces. The 
structure is subjected to prescribed displacements along the edges of the plates. 

The discretization and the decomposition into substructures for this test case are illustrated in 



Fig. 11 The total number of DOFs involved was 12 10 6 , distributed among 10,600 substructures. 
The number of macroscopic DOFs was 310 5 , which would have made the direct resolution on a 
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Figure 9: The LATIN convergence curves (error criterion vs. the number of iterations) for several 
numbers of macroitcrations 



standard computer very inefficient. 29 processors each with 4 Gigabytes memory were used for this 
calculation, which led to a super-coarse grid problem of dimension 150 (6 unknowns per floating super- 
substructure). 

Figure [T2| shows the damage map in the composite bolted joint after 70 time steps. The nonlinear 
calculation corresponding to each time step was carried out until the LATIN error criterion got below 
10~ 2 , which occurred after an average of 80 LATIN iterations per time step. The average CPU time 
required for the calculation of each time step was 30 minutes, which is reasonable considering the small 
number of processors used. 

However, the number of global LATIN iterations was quite high compared to what could have been 
obtained using the relocalization strategy of [TTj in the vicinity of the crack's front (see the results of 
Fig. [6]). Figure 12 is a clear illustration of the difficulties which may arise in the use of this dedicated 
technique in a general case in which multiple crack front propagations may be involved. The front has 
a complex shape, which raises the difficult issue of the choice of the number of relocalization zones and 
their sizes. In addition, this test case is very unstable. These instabilities were handled globally using 
an arc-length algorithm along with the three-scale resolution strategy, but local instabilities might 
also appear within the region extracted for the relocalization calculations, a situation which has not 
yet been addressed at this stage in our development. Therefore, in the future, it might be necessary 
to generalize the relocalization strategy in order to improve the efficiency of the enhanced multiscale 
domain decomposition technique for complex laminated structures. 



6 Conclusion 

The accurate prediction of delamination in extended process zones of laminated composite structures 
requires refined models of the material's behavior, leading to the resolution of huge systems of equa- 
tions. In order to solve such problems accurately, we used a two-scale domain decomposition strategy 
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Figure 10: Composite bolted joint made of 16 0.125 mm-thick plies. The in-plane dimensions of the 
structure are 30 x 5 mm. Prescribed displacements are applied along the left-hand side of the [0 9O2 0] s 
composite plate and along the right-hand sides of the [0 90] s composite plates. 




Figure 11: Discretization of the composite bolted joint (12 10 6 DOFs), decomposition into substruc- 
tures (10, 600 substructures) and assignment to processors (29 CPUs) 




Figure 12: Damage map of the cohesive interfaces of the composite bolted joint after the 70 time 
step of the quasi-static incremental analysis procedure 
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based on an iterative resolution algorithm. This method is particularly appropriate for laminated 
mesomodels, in which 3D and 2D entities are introduced separately. 

This strategy was improved in order to enable it to handle very large delamination problems. 
A systematic analysis of the features of the method on the different scales was performed. First, we 
showed that in the high-gradient zones the classical scale separation was insufficient to ensure numerical 
scalability. Therefore, we developed a subresolution procedure which preserves the numerical scalability 
of the crack propagation analysis, but still needs to be automated for complex structures and regulated 
against local instabilities. We also proved that a third scale is required. Then, the problem on the 
intermediate scale was solved using a parallel iterative algorithm which enabled the rapid transmission 
of the very-large-wavelength part of the solution. Global instabilities were handled through a classical 
arc-length algorithm with local control (e.g. based on the maximum damage increment) and adjustment 
of the "time" steps during the calculation of the evolution of damage. 

In future developments, 3D analysis in the process zone will be used in conjunction with plate 
analysis, which would be sufficient to describe the solution in the low-gradient zones. We will also, 
within the MAAXIMUS project, investigate the interaction between delamination and buckling for the 
simulation of components of aeronautical structures. 
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